{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 3\n",
    "\n",
    "### <center> Автор материала: Дмитрий Коргун"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Индивидуальный проект по анализу данных поездок на велосипедах."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**План исследования**\n",
    " - Описание набора данных и признаков\n",
    " - Первичный анализ признаков\n",
    " - Первичный визуальный анализ признаков\n",
    " - Закономерности, \"инсайты\", особенности данных\n",
    " - Предобработка данных\n",
    " - Создание новых признаков и описание этого процесса\n",
    " - Кросс-валидация, подбор параметров\n",
    " - Построение кривых валидации и обучения \n",
    " - Прогноз для тестовой или отложенной выборки\n",
    " - Оценка модели с описанием выбранной метрики\n",
    " - Выводы\n",
    " \n",
    " Подробности на [kaggle](https://www.kaggle.com/benhamner/sf-bay-area-bike-share).</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime as datetime\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.model_selection import train_test_split, TimeSeriesSplit, cross_val_score, GridSearchCV\n",
    "from sklearn.model_selection import validation_curve, learning_curve\n",
    "from sklearn.metrics import r2_score, mean_absolute_error\n",
    "from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\n",
    "from sklearn.linear_model import LinearRegression, Lasso, Ridge\n",
    "\n",
    "from xgboost import XGBRegressor\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "from plotly import __version__\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "from plotly import graph_objs as go\n",
    "init_notebook_mode(connected = True)\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Интересующие нас данные хранятся в следующих файлах:\n",
    "1. trip.csv - данные об использовании великов\n",
    "    * **id** - идентификатор\n",
    "    * **duration** - продолжительность использования(секунд)\n",
    "    * **start_date** - дата начала использования велика\n",
    "    * **start_station_name** - имя пункта отправления\n",
    "    * **start_station_id** - идентификатор пункта отправления\n",
    "    * **end_date** - дата возвращения велика\n",
    "    * **end_station_name** - имя конечного пункта\n",
    "    * **end_station_id** - идентификатор конечного пункта\n",
    "    * **bike_id** - идентификатор велика\n",
    "2. weather.csv - данные о погодных условиях \n",
    "    * **date** - дата\n",
    "    * ***_max_temperature_f** - максимальная, средняя и минимальные температуры в °F\n",
    "    * ***_dew_point_f** - максимальная, средняя и минимальная температура точки росы\n",
    "    * ***_humidity** - максимальная, средняя и минимальная влажность\n",
    "    * ***_sea_level_pressure_inches** - максимальное, среднее и минимальное атмосферное давление\n",
    "    * ***_visibility_miles** - максимальная, средняя и минимальная видимость\n",
    "    * ***_wind_speed_mph** - максимальная и средняя скорость ветра\n",
    "    * **max_gust_speed_mph** - наибольшая скорость порывов ветра\n",
    "    * **precipitation_inches** - атмосферные осадки\n",
    "    * **cloud_cover** - облачность\n",
    "    * **events** - погодные условия (туман, дождь и т.д.)\n",
    "    * **wind_dir_degrees** - угол и направления ветра\n",
    "3. station.csv - информация о велосипедных парковках\n",
    "    * **id** - идентификатор\n",
    "    * **name** - имя\n",
    "    * **lat** - широта\n",
    "    * **long** - долгота\n",
    "    * **dock_count** - количество мест на парковке\n",
    "    * **city** - город\n",
    "    * **installation_date** - дата открытия\n",
    "    \n",
    "Целевой признак - количество поездок на великах в день"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 2. Первичный анализ признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загружаем данные из файлов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip = pd.read_csv('trip.csv', parse_dates=['start_date', 'end_date'])\n",
    "weather = pd.read_csv('weather.csv', parse_dates=['date'])\n",
    "station = pd.read_csv('station.csv', parse_dates=['installation_date'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Переведём продолжительность в минуты и посмотрим на распределение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip['duration'] //= 60\n",
    "trip[['duration']].describe(percentiles=[.10, .25, .50, .75, .90])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ого... Максимальное значение почти 200 дней, хватит на пару кругостветок (если вы, например, [этот](https://lenta.ru/news/2017/09/19/velorider/) парень). Думаю, что использование дольше нескольких дней можно смело отнести к выбросам, т.к. нас будут интересовать постоянные пользователи велопарковки.\n",
    "\n",
    "Как правило люди берут велик, что бы доехать от одной парковки до другой, т.к. 90% пользуются им меньше 20 минут, чего (по моему сугубо личному мнению) маловато, чтобы насладиться поездкой, т.е. велосипеды используются как общественный транспорт, соответвенно спрос должен коррелировать с началом и окончанием рабочего дня и падать по выходным, позже мы проверим эту теорию. \n",
    "\n",
    "Посмотрим как выглядят наши данные"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим пропуски"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "У некоторых записей не заполнен почтовый код, но это не страшно, он нам не пригодится.\n",
    "\n",
    "Дальше разберёмся с погодными условиями, посмотрим на распределения погодных условий, данных о ней много, даже слишком. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather[[col for col in weather.columns if 'min' in col]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather[[col for col in weather.columns if 'mean' in col]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather[[col for col in weather.columns if 'max' in col]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим пропуски"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С данными о погоде всё не так радужно, ноо... лучше, чем могло бы быть.\n",
    "Наибольшее количество пропусков в признаке events. Это можно объяснить отсутствием особенных погодных условий, то есть был обычный ясный день. \n",
    "\n",
    "Сравним даты данных о погодных условиях и поездках"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.date.min(), weather.date.max(), weather.date.max() - weather.date.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip.start_date.min(), trip.end_date.max(), trip.end_date.max() - trip.start_date.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Промежутки совпадают, ничего лишнего нет, отлично."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Напоследок данные о пунктах выдачи"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "station.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip.start_station_id.unique().shape[0] == station.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Всего имеем 70 парковок и из каждой из них имеется информация в данных. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "station.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 3. Первичный визуальный анализ признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим количество поездок за каждый день"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip_days = trip['start_date'].apply(lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()\n",
    "trip_days.columns = ['trips']\n",
    "\n",
    "def plotMovingAverage(df, n,  title):\n",
    "    \n",
    "    rolling_mean = df.rolling(window=n).mean()\n",
    "    \n",
    "    trace0 = go.Scatter(\n",
    "        x = df.index,\n",
    "        y = df['trips'],\n",
    "        mode = 'lines',\n",
    "        name='Actual values'\n",
    "    )\n",
    "    \n",
    "    trace1 = go.Scatter(\n",
    "        x = rolling_mean.index,\n",
    "        y = rolling_mean['trips'],\n",
    "        mode = 'lines',\n",
    "        name='Rolling mean trend'\n",
    "    )\n",
    "    \n",
    "    data = [trace0, trace1]\n",
    "    layout = {'title': title}\n",
    "    fig = dict(data=data, layout=layout)\n",
    "    iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(trip_days, 7, 'Сглаживаем по неделям') \n",
    "plotMovingAverage(trip_days, 30, 'Сглаживаем по месяцам')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно заметить сразу несколько интересных вещей: \n",
    "* Во-первых, что очевидно, поездки зависят от времени года, минимальное количество зимой, максимальное - летом.\n",
    "\n",
    "* Во-вторых, похоже наибольшей популярность велики пользуется в будние дни, и падает на выходных. Убедимся в этом построим гистограмму распределения по дням недели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weekdays = {\n",
    "    0: 'Понедельник',\n",
    "    1: 'Вторник',\n",
    "    2: 'Среда',\n",
    "    3: 'Четверг', \n",
    "    4: 'Пятница', \n",
    "    5: 'Суббота', \n",
    "    6: 'Воскресенье'\n",
    "}\n",
    "\n",
    "trip_weekdays = trip['start_date'].apply(lambda dt: dt.weekday()).value_counts().sort_index().to_frame()\n",
    "trip_weekdays.columns = ['trips']\n",
    "trip_weekdays.index = trip_weekdays.index.map(lambda day: weekdays[day])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = go.Bar(\n",
    "    x = trip_weekdays.index,\n",
    "    y = trip_weekdays['trips']\n",
    ")\n",
    "\n",
    "data = [trace]\n",
    "layout = {'title': 'Распределение по дням недели'}\n",
    "fig = dict(data=data, layout=layout)\n",
    "iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip['is_weekend'] = trip['start_date'].apply(lambda dt: int(dt.weekday() >= 5))\n",
    "\n",
    "trip_weekday = trip[trip['is_weekend'] == 0]['start_date'].apply(\n",
    "    lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()\n",
    "trip_weekday.columns = ['trips']\n",
    "trip_weekend = trip[trip['is_weekend'] == 1]['start_date'].apply(\n",
    "    lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()\n",
    "trip_weekend.columns = ['trips']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace0 = go.Scatter(\n",
    "    x = trip_weekday.index,\n",
    "    y = trip_weekday['trips'],\n",
    "    mode = 'markers',\n",
    "    name='Будни'\n",
    ")\n",
    "\n",
    "trace1 = go.Scatter(\n",
    "    x = trip_weekend.index,\n",
    "    y = trip_weekend['trips'],\n",
    "    mode = 'markers',\n",
    "    name='Выходные'\n",
    ")\n",
    "\n",
    "data = [trace0, trace1]\n",
    "layout = {'title': 'Распределение по будням/выходным'}\n",
    "fig = dict(data=data, layout=layout)\n",
    "iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь рассмотрим распределение поездок по часам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip_hours = trip['start_date'].apply(lambda dt: dt.hour).value_counts().sort_index().to_frame()\n",
    "trip_hours.columns = ['trips']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = go.Bar(\n",
    "    x = trip_hours.index,\n",
    "    y = trip_hours['trips']\n",
    ")\n",
    "\n",
    "data = [trace]\n",
    "layout = {'title': 'Распределение по часам'}\n",
    "fig = dict(data=data, layout=layout)\n",
    "iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 4. Закономерности, \"инсайты\", особенности данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Итак по результатам предварительного анализа можно сделать следующие выводы:\n",
    "* Пропусков в данных нет (кроме почтового индекса)\n",
    "* Можно избавиться от выбросов, когда время поездки более суток\n",
    "* Есть недельный и сезонный тренды\n",
    "* Чаще люди пользуются великами в будни, чтобы добраться до работы или с работы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 5. Предобработка данных "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для начала посмотрим на количество выбросов и отфильтруем поездки более суток. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"Выбросы: {0:d} шт, {1:.2f}%\".format((trip['duration'] >= 60 * 24).sum(),\n",
    "                                     (trip['duration'] >= 60 * 24).sum() / trip.shape[0] * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trip = trip[trip['duration'] < 60 * 24]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь начнём формирование обучающей выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Просуммируем количество поездок по дням и отсортируем по дате\n",
    "data = trip['start_date'].apply(lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()\n",
    "\n",
    "# Приведём наименование\n",
    "data.columns = ['target']\n",
    "data['date'] = data.index\n",
    "data.reset_index(drop=True, inplace=True)\n",
    "#data = data.reindex(['date', 'target'], axis=1)\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь добавим информацию о погодных условиях\n",
    "\n",
    "Информация о погодных условиях разбита по почтовым кодам, так что очевидно данных о погоде больше, чем нужно"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape[0], weather.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather['zip_code'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на пропуски в погодных данных для разных кодов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for code in weather['zip_code'].unique():\n",
    "    print(\"Код: {0}. Количество пропусков: {1}\".format(code, weather[weather.zip_code == code].isnull().sum().values.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные с кодом 94107 имееют меньше всего пропусков. Велопарковки расположены относительно недалеко друг от друга, поэтому можно выбрать эти данные в качестве основных погодных условий без большого риска в потере качества предсказаний."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather = weather[weather['zip_code'] == 94107]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь разберёмся с пропусками в погодных данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как уже говорилось выше, пропуск значения признака events означает отсутсвие особенных погодных условий. Заполним все пропуски этого признака новым погодным условием - `Fair`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.loc[weather.events.isnull(), 'events'] = \"Fair\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather['events'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Среди значений погодных условий есть дублирующие `Rain` и `rain`. Оставим только первый вариант простой заменой одного на другое."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.loc[weather['events'] == 'rain', 'events'] = \"Rain\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отлично! <s>Сам себя не похвалишь...</s> Теперь уберём пропуски в признаке `Максимальная скорость порыва ветра`. \n",
    "\n",
    "Предположим, что скорость максимального порыва должна сильно коррелировать с максимальной скоростью ветра. Найдём медианные значения порыва ветра по скорости и заменим пропуски этим значением."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gust_by_wind = weather.groupby('max_wind_Speed_mph')['max_gust_speed_mph'].median()\n",
    "def fill_gust_by_wind(row):\n",
    "    row['max_gust_speed_mph'] = gust_by_wind[row['max_wind_Speed_mph']]\n",
    "    return row\n",
    "\n",
    "weather[weather['max_gust_speed_mph'].isnull()] = \\\n",
    "    weather[weather['max_gust_speed_mph'].isnull()].apply(fill_gust_by_wind, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на типы данных в колонках информации о погоде"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признак `precipitation_inches` имеет тип object, что немного противоречит природе вещей и крайне не нравится методам обучения, которые мы будем использовать. Исправим это и снова заполним получившиеся пропуски медианными значениями."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather['precipitation_inches'] = pd.to_numeric(weather['precipitation_inches'], errors = 'coerce')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather.loc[weather['precipitation_inches'].isnull(), 'precipitation_inches'] = \\\n",
    "    weather[weather['precipitation_inches'].notnull()]['precipitation_inches'].median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Остаётся только добавить информацию о велопарковках и собрать всё вместе в обучающую выборку."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "station.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посчитаем количество парковочных мест на каждый день из обучающей выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_docks = []\n",
    "for day in data['date']:\n",
    "    total_docks.append(sum(station[station.installation_date <= day].dock_count))\n",
    "data['total_docks'] = total_docks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 6. Создание новых признаков и описание этого процесса\n",
    "Теперь, когда мы закончили предварительную обработку данных сформируем признаки для обучающей выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weather = pd.concat([weather, pd.get_dummies(weather['events'])], axis=1)\n",
    "# Выкидываем лишние признаки и добавляем в обучающую выборку\n",
    "weather.drop(['events', 'zip_code'], axis=1, inplace=True)\n",
    "data = data.merge(weather, left_on='date', right_on='date')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так же добавим дополнительные временны́е признаки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['weekend'] = data['date'].apply(lambda dt: int(dt.weekday() >= 5))\n",
    "data['weekday'] = data['date'].apply(lambda dt: int(dt.weekday() < 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['year'] = data['date'].apply(lambda dt: dt.year % 2013)\n",
    "data['month'] = data['date'].apply(lambda dt: dt.month)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И наконец-таки сформируем выборку и ответы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = data.drop(['target', 'date'], axis=1), data['target']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_df = pd.concat([X, y], axis=1)\n",
    "full_df.rename(mapper={0: 'target'}, axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь, когда мы сформировали обучабщую выборку, давайте посмотрим на матрицу корреляций."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplots(figsize=(14,10))\n",
    "sns.heatmap(full_df.corr(), cmap=\"BuPu\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как и следовало ожидать, у погодных признаков одной категории высокая корреляция (напрмер, максимальная, средняя и минимальные температуры будут коррелировать). Это может спровоцировать линейные моделей на огромные весовые коэффициенты, но мы всегда можем избавиться от лишних признаков или воспользоваться регуляризацией.\n",
    "\n",
    "Так же можно заметить, что целевой коррелирует с признаками \"Будни\" и \"Выходные\", причем у признака \"Будни\" корреляция выше.\n",
    "Ещё имеем отрицательную корреляцию с признаком \"Rain-Thunderstorm\" (т.е. дождь с грозой), что тоже вполне логично.\n",
    "\n",
    "Соберём все данные и отложим последние три месяца для теста качества модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_size = X.shape[0] - 90\n",
    "X_train, y_train = X[:train_size], y[:train_size]\n",
    "X_test, y_test = X[train_size:], y[train_size:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 7. Кросс-валидация, подбор параметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Как всем известно, типичные метрики в задачах регрессии это средняя абсолютная (MAE) и среднеквадратичная (MSE) ошибки.\n",
    "MSE сильнее штрафует за большие отклонения и хорошо подходит для контроля качества во время обучения, но не позволяет сделать выводы о качестве полученного решения.\n",
    "\n",
    "Ещё есть `R2` мера:\n",
    "$$\\large \n",
    "R^2 = 1 - \\frac{\\sum_{i=1}^l (a(x_i) - y_i^2)}{\\sum_{i=1}^l (y_i - \\hat{y})^2}$$\n",
    "\n",
    "Чем ближе `R2` к единице, тем лучше модель объясняет данные.<br/>Если `R2` близка к нулю, предсказания близки к константным.<br>\n",
    "Фактически данная мера - нормированная среднеквадратичная ошибка. Можно сказать получаем интерпретируемую MSE. Её мы и будем использовать в данной задаче. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def score_model(clf, title):\n",
    "    tscv = TimeSeriesSplit(n_splits=10)\n",
    "    cv_res = np.median(cross_val_score(clf, X_train, y_train, cv=tscv, scoring='r2'))\n",
    "    print(title, \"R2: {0:.3f}\".format(cv_res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Попробуем обучить несколько линейных и \"деревянных\" моделей со стандартными параметрами(пусть это и не совсем честно), проведём кросс-валидацию по TimeSeriesSplit разбиению и проверим качество на отложенной выборке для того, что бы определиться с моделью.\n",
    "\n",
    "Использование TimeSeriesSplit позволяет получить последовательные обучающие выборки из прошлого и тестовую из будущего, что приближает задачу к реальной. Первое разбиение TimeSeriesSplit будет очень маленьким и качество обучения на нём, соответсвенно, будет негативно сказывается на среднем качестве по кросс-валидации, поэтому вместо среднего будем рассматривать медиану.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models = [\n",
    "    (LinearRegression(), \"Линейная регрессия\"),\n",
    "    (Ridge(random_state=17), \"Регрессия с регуляризатором L2\"),\n",
    "    (Lasso(random_state=17), \"Регрессия с регуляризатором L1\"),\n",
    "    (RandomForestRegressor(random_state=17), \"Случайный лес\"),\n",
    "    (GradientBoostingRegressor(random_state=17), \"Градиентный бустинг(sklearn)\"),\n",
    "    (XGBRegressor(random_state=17), \"Градиентный бустинг(xgb)\")\n",
    "]\n",
    "\n",
    "for pair in models:\n",
    "    score_model(*pair)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель линейной регрессии показала самый слабый результат. Это можно объяснить как следствие мультиколлинеарности.\n",
    "Регрессия с регуляризацией, показала схожий результат. Lasso и Ridge могут лучше, если подобрать параметр регуляризации.\n",
    "Случайный лес и бустинги имеют сравнимые результаты.\n",
    "\n",
    "Выберем в качестве основной модели xgb (у него всё таки лучший скор на кросс валидации) и подберём параметры."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb = XGBRegressor(random_state=17)\n",
    "tscv = TimeSeriesSplit(n_splits=10)\n",
    "parameters = {'n_estimators': range(10, 311, 30),\n",
    "              'max_depth': [3, 9, 15, 45],\n",
    "             }\n",
    "grid = GridSearchCV(xgb, parameters, n_jobs=5, \n",
    "                    cv=tscv, scoring='r2', refit=True)\n",
    "grid.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 8. Построение кривых валидации и обучения "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb = XGBRegressor(random_state=17, **grid.best_params_)\n",
    "param_range = list(range(10, 200, 10))\n",
    "tscv = TimeSeriesSplit(n_splits=10)\n",
    "train_scores, test_scores = validation_curve(\n",
    "    xgb, X_train, y_train,\n",
    "    param_name=\"n_estimators\", param_range=param_range,\n",
    "    cv=tscv, scoring='r2', n_jobs=-1)\n",
    "train_scores_mean = np.mean(train_scores, axis=1)\n",
    "train_scores_std = np.std(train_scores, axis=1)\n",
    "test_scores_mean = np.mean(test_scores, axis=1)\n",
    "test_scores_std = np.std(test_scores, axis=1)\n",
    "\n",
    "layout = go.Layout(title=\"Кривые валидации\",\n",
    "                      xaxis=dict(title=\"n_estimators\"),# type='log'),\n",
    "                      yaxis=dict(title=\"Score\"))\n",
    "\n",
    "lw = 2\n",
    "p1 = go.Scatter(x=param_range, y=train_scores_mean,\n",
    "                name=\"Training score\",\n",
    "                mode='lines', \n",
    "                line=dict(color=\"orange\", width=lw))\n",
    "\n",
    "p2 = go.Scatter(x=param_range, y=train_scores_mean - train_scores_std,\n",
    "                mode='lines', showlegend=False,\n",
    "                line=dict(color=\"orange\", width=1))\n",
    "\n",
    "p3 = go.Scatter(x=param_range, y=train_scores_mean + train_scores_std,\n",
    "                mode='lines', showlegend=False,\n",
    "                line=dict(color=\"orange\", width=1),\n",
    "                fill='tonexty')\n",
    "\n",
    "p4 = go.Scatter(x=param_range, y=test_scores_mean,\n",
    "                name=\"Cross-validation score\",\n",
    "                mode='lines', \n",
    "                line=dict(color=\"navy\", width=lw))\n",
    "\n",
    "p5 = go.Scatter(x=param_range, y=test_scores_mean - test_scores_std,\n",
    "                mode='lines', showlegend=False,\n",
    "                line=dict(color=\"navy\", width=1)) \n",
    "\n",
    "p6 = go.Scatter(x=param_range, y=test_scores_mean + test_scores_std,\n",
    "                mode='lines', showlegend=False,\n",
    "                line=dict(color=\"navy\", width=1),\n",
    "                fill='tonexty') \n",
    "\n",
    "fig = go.Figure(data=[p2, p3, p5, p6, p1, p4], layout=layout)\n",
    "iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "После 50 деревьев увеличение сложности по данному параметру не улучшает качество модели, но кривые валидации всё ещё находятся достаточно далеко друг от друга. Следовательно, увеличение сложности модели, например, по другим параметрам, вероятно повысит качество модели. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_learning_curve(estimator, X, y, cv=None,\n",
    "                        n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):\n",
    "    train_sizes, train_scores, test_scores = learning_curve(\n",
    "        estimator, X, y, \n",
    "        cv=cv, n_jobs=n_jobs, train_sizes=train_sizes,\n",
    "        scoring='r2', random_state=17\n",
    "    )\n",
    "    train_scores_mean = np.mean(train_scores, axis=1)\n",
    "    train_scores_std = np.std(train_scores, axis=1)\n",
    "    test_scores_mean = np.mean(test_scores, axis=1)\n",
    "    test_scores_std = np.std(test_scores, axis=1)\n",
    "    \n",
    "    layout = go.Layout(title=\"Кривые обучения\")\n",
    "    \n",
    "    p1 = go.Scatter(x=train_sizes, y=test_scores_mean + test_scores_std,\n",
    "                    mode='lines',\n",
    "                    line=dict(color=\"green\", width=1),\n",
    "                    showlegend=False)\n",
    "    \n",
    "    p2 = go.Scatter(x=train_sizes, y=test_scores_mean - test_scores_std,\n",
    "                    mode='lines',\n",
    "                    line=dict(color=\"green\", width=1),\n",
    "                    showlegend=False, fill='tonexty')\n",
    "    \n",
    "    p3 = go.Scatter(x=train_sizes, y=train_scores_mean + train_scores_std,\n",
    "                    mode='lines',\n",
    "                    line=dict(color=\"red\", width=1),\n",
    "                    showlegend=False)\n",
    "    \n",
    "    p4 = go.Scatter(x=train_sizes, y=train_scores_mean - train_scores_std,\n",
    "                    mode='lines',\n",
    "                    line=dict(color=\"red\", width=1),\n",
    "                    showlegend=False, fill='tonexty')\n",
    "\n",
    "    p5 = go.Scatter(x=train_sizes, y=train_scores_mean, \n",
    "                    marker=dict(color='red'),\n",
    "                    name=\"Training score\", showlegend=True)\n",
    "    \n",
    "    p6 = go.Scatter(x=train_sizes, y=test_scores_mean, \n",
    "                    marker=dict(color='green'),\n",
    "                    name=\"Cross-validation score\", showlegend=True)\n",
    "    \n",
    "    fig = go.Figure(data=[p1, p2, p3, p4, p5, p6], layout=layout)\n",
    "    iplot(fig, show_link=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb = XGBRegressor(random_state=17, **grid.best_params_)\n",
    "tscv = TimeSeriesSplit(n_splits=10)\n",
    "plot_learning_curve(xgb, X_train, y_train, cv=tscv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кривые обучения не сошлись, информации недостаточно и добавление новых данных может улучшить результат."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 9. Прогноз для тестовой или отложенной выборки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При формировании обучающей выборки, мы сделали отложенную выборку из трёх последних месяцев, проверим насколько хороша наша модель."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb = XGBRegressor(random_state=17, **grid.best_params_)\n",
    "xgb.fit(X_train, y_train)\n",
    "y_pred = xgb.predict(X_test)\n",
    "print(\"R2: {0:.3f}, MAE: {1:.3f}\".format(r2_score(y_test, y_pred),\n",
    "                                         mean_absolute_error(y_test, y_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 10. Оценка модели с описанием выбранной метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Мы получили достаточно неплохие результаты R2 - 0.84, т.е. модель хорошо объясняет данные, MAE - 108 т.е. модель ошибается в среднем всего на 100 велосипедов, тоже приемлимо учитывая масштабы целевой переменной. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.median(), y_test.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так же можно посмотреть на важность признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importances = xgb.feature_importances_\n",
    "indices = np.argsort(importances)\n",
    "\n",
    "plt.subplots(figsize=(14,10))\n",
    "plt.title('Важность признаков')\n",
    "plt.barh(range(len(indices)), importances[indices], color='b', align='center')\n",
    "plt.yticks(range(len(indices)), X_train.columns[indices])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 11. Выводы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Получили интересную модель предсказывающую количество поездок на велосипедах по погодным условиям. Самые важные, по мнению модели, признаки это:\n",
    "* Температура\n",
    "* Ветер\n",
    "* Количество осадков\n",
    "Качество модели на тестовой выборке по метрике R2 составляет 0.838, что означает хорошую обобщающую способность. \n",
    "\n",
    "Возможные улучшения:\n",
    "* Как показали кривые обучения, добавление новых должно повысить качество модели\n",
    "* Избавиться от мультиколлинеарнсоти признаков (регуляризацией, PCA или отбором признаков)\n",
    "* Создание новых признаков или использование открытых данных\n",
    "* Более тонкая настройка параметров"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
